#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <cusolverDn.h>
#include <stdio.h>
#include"library1.h"
#include<string.h>
#include<stdlib.h>
#include <time.h>
#include<math.h>      //gcc�����ʱ����Ҫ��-lmѡ��

#define THREASHOLD 1e-8
#define ITERATION 11
#define COL 1
struct node
{
	int min;
	int max;
};
__device__ int sign(double number) {
    if(number<0)
        return -1;
    else
        return 1;
}

int procSize;   //���������Ŀ

/*���ļ��ж�ȡ����*/
void readFromFile(double *matrix,int row,int col,char* file){
    FILE *fp;
    fp=fopen(file,"r");//���ļ�
    for(int i=0;i<row;i++) //������
        for(int j=0;j<col;j++)
            fscanf(fp,"%lf",&matrix[i*col+j]);
    fclose(fp);//�ر��ļ�
}

/*�Ѿ���д���ļ�*/
void writeToFile(double *matrix,int rows,int columns,char* file){
    FILE *fp;
    if((fp=fopen(file,"w"))==NULL){
        perror("fopen");
        exit(1);
    }
    fprintf(fp,"%d\t%d\n",rows,columns);        //�ļ����м�¼���������������
    int i,j;
    for(i=0;i<rows;++i){
        for(j=0;j<columns;++j){
            fprintf(fp,"%-10f\t",matrix[i*columns+j]);
        }
        fprintf(fp,"\n");
    }
    fclose(fp);
}

/*������д���ļ�*/
void vectorToFile(double *vector,int len,char* file){
    FILE *fp;
    if((fp=fopen(file,"w"))==NULL){
        perror("fopen");
        exit(1);
    }
    int i;
    for(i=0;i<len;++i){
        fprintf(fp,"%-10f\t",vector[i]);
    }
    fclose(fp);
}

/*�����������е���Jacobi�����任*/
__device__ void orthogonalVector(double *Ci,double *Cj,int len1,int i,int j,double *d,double *Vi,double *Vj,int len2,int *pass){
    double ele=vectorProduct(Ci,Cj,len1);
	printf("threadIdx=%d,ele=%lf\n",threadIdx.x,ele);
    if(fabs(ele)<THREASHOLD)
        return;     //��������Ѿ�����������Ҫ���б任���򷵻�true
    *pass=0;
    double ele1=d[i];
    double ele2=d[j];


    double tao=(ele1-ele2)/(2*ele);
	//printf("threadIdx=%d,tao=%lf\n",threadIdx.x,tao);
    double tan=sign(tao)/(fabs(tao)+sqrt(1+pow(tao,2)));
	//printf("threadIdx=%d,tan=%lf\n",threadIdx.x,tan);
    double cos=1/sqrt(1+pow(tan,2));
	//printf("threadIdx=%d,cos=%lf\n",threadIdx.x,cos);
    double sin=cos*tan;
	//printf("threadIdx=%d,sin=%lf\n",threadIdx.x,sin);
    int row;
    for(row=0;row<len1;++row){
		//printf("threadIdx=%d,Ci=%lf\n",threadIdx.x,Ci[row]);
		//printf("threadIdx=%d,Cj=%lf\n",threadIdx.x,Cj[row]);
        double var1=Ci[row]*cos+Cj[row]*sin;
		//printf("threadIdx=%d,var1=%lf\n",threadIdx.x,var1);
        double var2=Cj[row]*cos-Ci[row]*sin;
		//printf("threadIdx=%d,var2=%lf\n",threadIdx.x,var2);

        Ci[row]=var1;
        Cj[row]=var2;

    }
	d[i]=vectorProduct(Ci,Ci,len1);
	d[j]=vectorProduct(Cj,Cj,len1);
	//printf("d[i]=%lf\n",d[i]);
	//printf("d[j]=%lf\n",d[j]);
    for(row=0;row<len2;++row){
        double var1=Vi[row]*cos+Vj[row]*sin;
        double var2=Vj[row]*cos-Vi[row]*sin;

        Vi[row]=var1;
        Vj[row]=var2;

    }
}

/*__device__ void sort(double a[],int p[],int max)
{
	int min = 0;
	double key = a[min];
	int i = min;
	int j = max;
	double temp;
	int temp1;
    struct node Stack[100];
    int top = 0;
    Stack[top].min = min;
    Stack[top].max = max;

	while(top>-1)
	{
		//min max��¼��ǰ����������������޺��м���
		i = min = Stack[top].min;
		j = max = Stack[top].max;
		top--;
		key = a[min];

		while(i<j)
		{
			while((i<j) && (key <= a[j]))
			{j--;}
			if(i < j)
			{
				temp = a[i];a[i] = a[j];a[j] = temp;
				temp1 = p[i];p[i] = p[j];p[j] = temp1;
				i++;
			}

			while((i<j) && (key >= a[i]))
			{i++;}
			if(i < j)
			{
				temp = a[i];a[i] = a[j];a[j] = temp;
				temp1 = p[i];p[i] = p[j];p[j] = temp1;
				j--;
			} 
		}//����һ��  �����Ȱ�ֵС��ȫ�������  �Ȱ�ֵ��ķ��ұ�

		if(min < i-1)
		{
			top++;
			Stack[top].min = min;
			Stack[top].max = i-1;
		}
		if(max>i+1)
		{
			top++;
			Stack[top].min = i+1;
			Stack[top].max = max;
			}
	} 
}*/

__device__ void sort(double *a,int *p,int max){
	__shared__ double tnorm[700];
	__shared__ int tp[700];
	int tid=threadIdx.x;
	/*while(tid<max){
		printf("a[%d]=%lf,p[%d]=%d\n",tid,a[tid],tid,p[tid]);
		tid=tid+blockDim.x;
	}
	tid=threadIdx.x;*/
	int i,j,k;
	int n=1;  //�鲢��������Ԫ�ظ���
	while(n<max){
		while(tid<max){
			if((tid%(n*2)==0)){
				for(i=tid,j=tid+n,k=tid;;k++){
					if(k>=max)
						break;
					//printf("a[%d]=%lf\n",tid,a[tid]);
					if(a[i]>a[j]||j>=max){
						/*if(p[i]>500)
							printf("1:tid=%d,a[%d]=%lf,p[%d]=%d,tp[%d]=%d,%d\n",tid,i,a[i],i,p[i],k,tp[k],n);*/
						tnorm[k]=a[i];
						tp[k]=p[i];
						i++;
						if(i==tid+n){
							for(;j<tid+n+n;j++){
								k++;
								if(k>=max)
									break;
								/*if(p[i]>500)
									printf("2:tid=%d,a[%d]=%lf,p[%d]=%d,tp[%d]=%d,%d\n",tid,j,a[j],j,p[j],k,tp[k],n);*/
								tnorm[k]=a[j];
								tp[k]=p[j];
								
							}
							break;
						}
					}else{
						/*if(p[i]>500)
							printf("23:tid=,%d,a[%d]=%lf,p[%d]=%d,tp[%d]=%d,%d\n",tid,j,a[j],j,p[j],k,tp[k],n);*/
						tnorm[k]=a[j];
						tp[k]=p[j];
						j++;
						if(j==tid+n+n){
							for(;i<tid+n;i++){
								k++;
								if(k>=max)
									break;
								/*if(p[i]>500)
									printf("15:tid=%d,a[%d]=%lf,p[%d]=%d,tp[%d]=%d,%d\n",tid,i,a[i],i,p[i],k,tp[k],n);*/
								tnorm[k]=a[i];
								tp[k]=p[i];
							}
							break;
						}
					}
				}
			}
			tid=tid+blockDim.x;
		}
		n*=2;
		__syncthreads();
		tid=threadIdx.x;
		while(tid<max){
			/*if(tp[tid]>500){
				printf("a[%d]=%lf,p[%d]=%d,tp[%d]=%d,%d\n",tid,a[tid],tid,p[tid],tid,tp[tid],threadIdx.x);
			}*/
			a[tid]=tnorm[tid];
			p[tid]=tp[tid];
			/*if(p[tid]>500){
				printf("a[%d]=%lf,p[%d]=%d,tp[%d]=%d,%d,%d\n",tid,a[tid],tid,p[tid],tid,tp[tid],n,threadIdx.x);
			}*/
			tid=tid+blockDim.x;
		}
		tid=threadIdx.x;
		__syncthreads();
	}
}

__device__ void push_forward(int *p_a,int *p_b,int i,int length_b){
	volatile int temp;
	__shared__ int p[500];
	int tid=threadIdx.x;
	if(threadIdx.x==0){
		temp=p_a[length_b-i-1];
		p_a[length_b-i-1]=p_b[length_b-i-1];
		p_b[length_b-i-1]=temp;
		temp = p_b[0];
	}
	__syncthreads();
	while(tid<length_b-1){
		p[tid]=p_b[tid+1];
		tid+=blockDim.x;
	}
	__syncthreads();
	if(threadIdx.x==0)
		p[length_b-1]=temp;
	tid=threadIdx.x;
	__syncthreads();
	while(tid<length_b){
		p_b[tid]=p[tid];
		tid+=blockDim.x;
	}
	__syncthreads();
}

__device__ void push_forward1(int *p_a,int *p_b,int length_b){
	volatile int temp;
	__shared__ int p[500];
	int tid=threadIdx.x;
	if(threadIdx.x==0){
		temp=p_a[length_b];
	}
	__syncthreads();
	while(tid<length_b){
		p[tid+1]=p_a[tid];
		tid+=blockDim.x;
	}
	tid=threadIdx.x;
	__syncthreads();
	if(threadIdx.x==0)
		p[0]=p_b[0];
	__syncthreads();
	while(tid<=length_b){
		p_a[tid]=p[tid];
		tid+=blockDim.x;
	}
	tid=threadIdx.x;
	__syncthreads();
	while(tid<length_b-1){
		p[tid]=p_b[tid+1];
		tid+=blockDim.x;
	}
	tid=threadIdx.x;
	__syncthreads();
	if(threadIdx.x==0)
		p[length_b-1]=temp;
	__syncthreads();
	while(tid<length_b){
		p_b[tid]=p[tid];
		tid+=blockDim.x;
	}
	__syncthreads();
	}

__device__ int p_a[2][700];
__device__ int p1[700];
__device__ int p_ab[2][700][700];
__device__ int groupColumnIndex[500][50];  //[2*gridDim.x][col/(2*gridDim.x)];

__device__ void initilize_index_of_one_side_jacobi(int *p_a, int *p_b, int *i, int k,int length_b,int row){
	if(row%2==0){
		if(((k+1)==2*(*i+1)+1)&&threadIdx.x==0){
			*i=*i+1;
		}
		__syncthreads();
		push_forward(p_a,p_b,*i,length_b);
	}
	else{
		push_forward1(p_a,p_b,length_b);
	}

    /*
	if(threadIdx.x==0 && k < 498){
		
		for(int i=0;i<length_b;i++){
		
			if ( p_a[i] != p1[p_ab[0][i][k+1]] || p_b[i] != p1[p_ab[1][i][k+1]] )
			{
				printf("p_a[0][%d]=%d, p_ab[0][%d][%d]=%d\n", i, p_a[i],i, k+1, p1[p_ab[0][i][k+1]]);
				printf("p_a[1][%d]=%d, p_ab[1][%d][%d]=%d\n", i, p_b[i],i, k+1, p1[p_ab[1][i][k+1]]);
			}
		}
		
	}
    */

	__syncthreads();
}	

//���������ĳ����
__device__ void exchangeRow(double *matrix,int rows,int columns,int i,int j){
    assert(i<columns);
    int col;
    for(col=0;col<columns;++col){
        double tmp=matrix[i*columns+col];
        matrix[i*columns+col]=matrix[j*columns+col];
        matrix[j*columns+col]=tmp;
    }
}

//�����������ڻ�
__device__ void vectorProduct1(double *a,double *b,int len,double *cache){
    //ÿһ�����϶���cache�����Ŀ������໥֮�䲻Ӱ��
	cache[threadIdx.x]=0.0;
    //tidΪ�̵߳�ƫ����
    int tid=threadIdx.x;
    double temp=0.0;
	while(tid<len){
		temp+=a[tid]*b[tid];
		tid+=blockDim.x;
	}
	cache[threadIdx.x]=temp;
    //ͬ������ǰ���ϵ��߳�
    __syncthreads();

    //�ڿ��ڼ��㲿�ֺ�
	if(threadIdx.x<128)
		cache[threadIdx.x]+=cache[threadIdx.x+128];
	__syncthreads();
	if(threadIdx.x<64)
		cache[threadIdx.x]+=cache[threadIdx.x+64];
	__syncthreads();
	if(threadIdx.x<32){
		volatile double* mysh = &cache[threadIdx.x];
        *mysh += mysh[32];
        *mysh += mysh[16];
        *mysh += mysh[8];
        *mysh += mysh[4];
        *mysh += mysh[2];
        *mysh += mysh[1];
	}
	
	
	/*int lane=threadIdx.x&31;
	double lane_local = cache[threadIdx.x];
	double lane_recv;
    lane_recv = __shfl_xor_sync(0xffffffff,lane_local, 1);
	if ((lane & 1) == 0)
		lane_local += lane_recv;
    lane_recv = __shfl_xor_sync(0xffffffff,lane_local, 2);
	if ((lane & 3) == 0)
		lane_local += lane_recv;
    lane_recv = __shfl_xor_sync(0xffffffff,lane_local, 4);
	if ((lane & 7) == 0)
		lane_local += lane_recv;
    lane_recv = __shfl_xor_sync(0xffffffff,lane_local, 8);
	if ((lane & 0xF) == 0)
		lane_local += lane_recv;
    lane_recv = __shfl_xor_sync(0xffffffff,lane_local, 0x10);
	if (lane == 0){
		lane_local += lane_recv;
	}
	cache[threadIdx.x]=__shfl_sync(0xffffffff,lane_local,0);
	__syncthreads();*/
	
	/*for (int offset = 16; offset > 0; offset /= 2)
        lane_local += __shfl_down_sync(0xffffffff, lane_local, offset);
	cache[threadIdx.x]=__shfl_sync(0xffffffff,lane_local,0);
	__syncthreads();*/
	
	/*volatile double* mysh = &cache[threadIdx.x];                  //warp1;
	*mysh += mysh[16];
	*mysh += mysh[8];
	*mysh += mysh[4];
	*mysh += mysh[2];
	*mysh += mysh[1];
	__syncthreads();*/
	
	/*if (lane == 0){                   //wrap��Լ2
		cache[threadIdx.x>>5] = lane_local;
	}
	__syncthreads();
	
	lane_local = cache[lane & 7];
	lane_recv = __shfl_xor_sync(0xffffffff,lane_local, 1, 8);
	if ((lane & 1) == 0){
		lane_local += lane_recv;
	}
	lane_recv = __shfl_xor_sync(0xffffffff,lane_local, 2, 8);
	if ((lane & 3) == 0){
		lane_local += lane_recv;
	}
	lane_recv = __shfl_xor_sync(0xffffffff,lane_local, 4, 8);
	if (lane == 0){
		cache[threadIdx.x] = lane_local + lane_recv;
	}*/
	
	/*if(threadIdx.x<4){             //warp��Լ1 
		volatile double* mysh1 = &cache[threadIdx.x*32];
        *mysh1 += mysh1[4*32];
        *mysh1 += mysh1[2*32];
        *mysh1 += mysh1[32];
	}*/
}
__device__ void vectorProduct2(double *a,double *b,int len,double *cache,double *cache1){
    //ÿһ�����϶���cache�����Ŀ������໥֮�䲻Ӱ��
	cache[threadIdx.x]=0.0;
	cache1[threadIdx.x]=0.0;
    //tidΪ�̵߳�ƫ����
    int tid=threadIdx.x;
	while(tid<len){
		cache[tid]+=a[tid]*a[tid];
		//cache1[tid]+=b[tid]*b[tid];
		tid+=blockDim.x;
	}
    //ͬ������ǰ���ϵ��߳�
    __syncthreads();

    //�ڿ��ڼ��㲿�ֺ�
	if(threadIdx.x<128)
		cache[threadIdx.x]+=cache[threadIdx.x+128];
		//cache1[threadIdx.x]+=cache1[threadIdx.x+128];
	__syncthreads();
	if(threadIdx.x<64)
		cache[threadIdx.x]+=cache[threadIdx.x+64];
		//cache1[threadIdx.x]+=cache1[threadIdx.x+64];
	__syncthreads();
	if(threadIdx.x<32){
		volatile double* mysh = &cache[threadIdx.x];
		//volatile double* mysh1 = &cache1[threadIdx.x];
        *mysh += mysh[32];
		//*mysh1 += mysh1[32];
        *mysh += mysh[16];
		//*mysh1 += mysh1[16];
        *mysh += mysh[8];
		//*mysh1 += mysh1[8];
        *mysh += mysh[4];
		//*mysh1 += mysh1[4];
        *mysh += mysh[2];
		//*mysh1 += mysh1[2];
        *mysh += mysh[1];
		//*mysh1 += mysh1[1];
	}
}
__global__ void sort1(double *A,double *V,double *U,double *S,double *A1,double *V1,int row,int col,size_t pitch){
	__shared__ double d[700];
	__shared__ int p[700];
	int tid=threadIdx.x;
	while(tid<row){
		double* r=(double*)((char*)A+tid*pitch);
		d[tid]=sqrt(vectorProduct(r,r,row));
		p[tid]=tid;
		tid=tid+blockDim.x;
	}
	__syncthreads();
	sort(d,p,row);
	__syncthreads();
	tid=threadIdx.x;
	while(tid<row){
		S[tid]=d[tid];
		//printf("S[%d]=%lf,d[%d]=%lf,p[%d]=%d\n",threadIdx.x,S[threadIdx.x],threadIdx.x,d[threadIdx.x],threadIdx.x,p[threadIdx.x]);
		double* Ai=(double*)((char*)A1+tid*pitch);
		double* Aj=(double*)((char*)A+p[tid]*pitch);
		double* Vi=(double*)((char*)V1+tid*pitch);
		double* Vj=(double*)((char*)V+p[tid]*pitch);
		for(int i=0;i<row;i++){
			Ai[i]=Aj[i];
			Vi[i]=Vj[i];
		}
		tid=tid+blockDim.x;
	}
}
__global__ void getU(double *A1,double *S,double *U,int row,int col,int pitch){
	int tid=threadIdx.x;
	while(tid<row){
		if(S[tid]<THREASHOLD){
			for(int i=0;i<col;++i)
				U[tid*col+i]=0.0;
			}
		else
			for(int i=0;i<col;++i){
				double* Ai=(double*)((char*)A1+tid*pitch);
				U[tid*col+i]=Ai[i]/S[tid];
			}
		tid=tid+blockDim.x;
	}
}

__device__ volatile int Arrayin[500];
__device__ void __gpu_sync(volatile int *Arrayin)
{ 
	// only thread0 is used for synchonization
	if (threadIdx.x == 0)
	{
		Arrayin[blockIdx.x] = 1;
		__threadfence();
	}
 
	if (blockIdx.x == gridDim.x-1)
	{
		int tid=threadIdx.x;
		while (tid < gridDim.x)
		{
			while (Arrayin[tid] != 1)
			{
				// Do nothing here
			}
			tid=tid+blockDim.x;
		}
 
		__syncthreads();
		tid=threadIdx.x;
		while (tid < gridDim.x)
		{
			Arrayin[tid] = 2;
			tid=tid+blockDim.x;
		}
	}
 
	if (threadIdx.x == 0)
	{
		while (Arrayin[blockIdx.x] != 2)
		{
			// Do nothing here
		}
	}
 
	__syncthreads();
}

__device__ void __gpu_syncNew(volatile int *Arrayin, int iter)
{
	// only thread0 is used for synchronization
	if ( threadIdx.x == 0 )
	{

	    Arrayin[blockIdx.x] += 1;
		__threadfence();


		if ( blockIdx.x < gridDim.x-1 )
		{
			while ( Arrayin[blockIdx.x+1] < iter )
			{
				// Do nothing here
			}
		}
			
		if ( blockIdx.x == gridDim.x-1 )
		{
            while ( Arrayin[0] < iter )
            {
                // Do nothing here
            }
		}
	}

	__syncthreads();
}

__device__ double d1[700];
__device__ double d2[700];
__global__ void getD(double *A,int row,int col,size_t pitch){
	__shared__ double cache[512];
	double* r=(double*)((char*)A+blockIdx.x*pitch);
	vectorProduct1(r,r,row,cache);
	if(threadIdx.x==0){
		d1[blockIdx.x]=cache[0];
        d2[blockIdx.x]=cache[0];
	}
}
__global__ void getrank(int row,int col,int *pass){
	__shared__ int p[700];
	int tid=threadIdx.x;
	while(tid<row){
		pass[tid]=1;
		p[tid]=tid;
		tid=tid+blockDim.x;
	}
	__syncthreads();
	sort(d2,p,row);
	__syncthreads();
	if(threadIdx.x==0){
		int a=0,b=0;
		for(int i=0;i<row;i++){
			//printf("p[%d]=%d,d[%d]=%lf\n",i,p[i],i,d[i]);
			if(i%2==0){
				p_a[0][a]=p[i];
				//printf("p_a[0][%d]=%d, p_ab[0][%d][0]=%d\n", i, p_a[0][a],i, p_ab[0][a][0]);
				a++;
			}
			else{
				p_a[1][b]=p[i];
				//printf("p_a[1][%d]=%d, p_ab[1][%d][1]=%d\n", i, p_a[1][b],i, p_ab[1][b][0]);
				b++;
			}

			p1[i] = p[i];
		}
	}
	__syncthreads();
}

__global__ void getrankNew(int row,int col,int *pass){
	__shared__ int p[700];
	int tid=threadIdx.x;
	while(tid<row){
		pass[tid]=1;
		p[tid]=tid;;
		tid=tid+blockDim.x;
	}
	__syncthreads();
	sort(d1,p,row);
	__syncthreads();
	if(threadIdx.x==0){
		int h=0,k=0;
		for(int i=0;i<row;i++){
			if(i%2==0){
				p_ab[0][h][0]=p[i];
				p_a[0][h++]=p[i];
			}
			else{
				p_ab[1][k][0]=p[i];
				p_a[1][k++]=p[i];
			}
		}
	}
	__syncthreads();
	
	int count = 0;
	for ( int i = 1; i < row-(1-row%2); i ++ ){

		if(row%2==0){
			
			if( ( i-1 == 2*(count+1) ) && threadIdx.x == 0 ){
				count = count + 1;
			}

			__syncthreads();
			
			push_forward(p_a[0],p_a[1],count,row/2);
		}
		else{
			push_forward1(p_a[0],p_a[1],row/2);
		}
		
		__syncthreads();
		tid=threadIdx.x;
		while(tid<row/2){
			p_ab[0][tid][i]=p_a[0][tid];
			p_ab[1][tid][i]=p_a[1][tid];
			tid=tid+blockDim.x;
		}
		__syncthreads();
	
	}
	/*
	if (threadIdx.x == 0)
		for ( int k = 0; k < row-1; k ++ )
		{
			printf("k=%d\n",k);
			for ( int i = 0; i< 2; i ++ )
			{
				printf("i=%d\t,",i);
				for (int j = 0; j < row/2; j ++ )
				{
					printf("%d ", p_ab[i][j][k]);
				}
				printf("\n");
			}
			printf("\n");
		}
	__syncthreads();
	*/
}


__global__ void getRankNewNew(int row){

    int tid;

    if( threadIdx.x == 0 ){
        int h=0, k=0;
        for ( int i = 0; i < row; i ++ ){
            if ( i%2 == 0 ){
                p_ab[0][h][0]=i;
                p_a[0][h++]=i;
            }
            else{
                p_ab[1][k][0]=i;
                p_a[1][k++]=i;
            }
        }
    }
    __syncthreads();

    int count = 0;
    for ( int i = 1; i < row-(1-row%2); i ++ ) {

        if (row % 2 == 0) {

            if ((i - 1 == 2 * (count + 1)) && threadIdx.x == 0) {
                count = count + 1;
            }

            __syncthreads();

            push_forward(p_a[0], p_a[1], count, row / 2);
        }
        else {
            push_forward1(p_a[0], p_a[1], row / 2);
        }

        __syncthreads();
        tid = threadIdx.x;
        while (tid < row / 2) {
            p_ab[0][tid][i] = p_a[0][tid];
            p_ab[1][tid][i] = p_a[1][tid];

            tid = tid + blockDim.x;
        }
        __syncthreads();
    }
}

__global__ void gpusvd(double *A,double *V, int row, int col,size_t pitch,int *pass,int *count_i) {       //�Ǽ��ɲ���svd�ֽ�
    __shared__ double ele, ele1, ele2;
    __shared__ double cache[256];
    __shared__ double Acol[2][500];
    __shared__ double Vcol[2][500];
    __shared__ int p_index;
    int pa[2];
    int tid;

    if (threadIdx.x == 0){
        Arrayin[blockIdx.x] = 0;
    }
    __syncthreads();

	double *Aco[2],*Vco[2];
	for(int i=0;i<row-(1-row%2);i++){

	    //pa[0]=p_a[0][blockIdx.x];
		//pa[1]=p_a[1][blockIdx.x];

		pa[0] = p1[p_ab[0][blockIdx.x][i]];
        pa[1] = p1[p_ab[1][blockIdx.x][i]];

		// main computation

			Aco[0]=(double*)((char*)A+pa[0]*pitch);
			Aco[1]=(double*)((char*)A+pa[1]*pitch);
			Vco[0]=(double*)((char*)V+pa[0]*pitch);
			Vco[1]=(double*)((char*)V+pa[1]*pitch);

			if ( i == 0 ) {
                tid = threadIdx.x;
                while (tid < row) {
                    Acol[0][tid] = Aco[0][tid];
                    Acol[1][tid] = Aco[1][tid];
                    Vcol[0][tid] = Vco[0][tid];
                    Vcol[1][tid] = Vco[1][tid];
                    tid = tid + blockDim.x;
                }
                __syncthreads();
            }
			else {

                if(p_index==0) {
                    tid = threadIdx.x;
                    while(tid < row){
                        Acol[0][tid] = Acol[1][tid];
                        tid = tid + blockDim.x;
                    }
                    __syncthreads();
                }

                tid=threadIdx.x;
                while(tid<row){
                    Acol[1][tid]=Aco[1][tid];
                    Vcol[1][tid]=Vco[1][tid];
                    tid=tid+blockDim.x;
                }
				__syncthreads();
			}

            if(threadIdx.x==0){
                p_index=(blockIdx.x!=(row/2-1-i/2));
            }
            __syncthreads();

			vectorProduct1(Acol[0],Acol[1],row,cache); //5

			if(threadIdx.x==0){
				ele=cache[0];
			}
			__syncthreads();

		if (fabs(ele)>THREASHOLD) {
            if (threadIdx.x == 0) {
                ele1 = d1[pa[0]];
                ele2 = d1[pa[1]];
                pass[blockIdx.x] = 0;
            }
            __syncthreads();

            double aro;
            double aco;
            double tao = (ele1 - ele2) / (2 * ele);
            double tan = sign(tao) / (fabs(tao) + sqrt(1 + pow(tao, 2)));
            double cos = 1 / sqrt(1 + pow(tan, 2));
            double sin = cos * tan;

            /*
            int j;
            if ( p_index == 0 ){
                j = 0;
            }
            else{
                j = 1;
            }
            __syncthreads();

            tid = threadIdx.x;
            while (tid < row) {
                aro = Acol[0][tid] * cos + Acol[1][tid] * sin;
                aco = Acol[1][tid] * cos - Acol[0][tid] * sin;
                Acol[0][tid] = aro;
                Acol[1][tid] = aco;
                aro = Vcol[0][tid] * cos + Vcol[1][tid] * sin;
                aco = Vcol[1][tid] * cos - Vcol[0][tid] * sin;
                Vcol[0][tid] = aro;
                Vcol[1][tid] = aco;

                Aco[j][tid] = Acol[j][tid];
                Vco[j][tid] = Vcol[j][tid];
                tid = tid + blockDim.x;
            }
            __syncthreads();
            */

            if ( p_index == 0 ){
                tid = threadIdx.x;
                while (tid < row) {
                    aro = Acol[0][tid] * cos + Acol[1][tid] * sin;
                    aco = Acol[1][tid] * cos - Acol[0][tid] * sin;
                    Acol[0][tid] = aro;
                    Acol[1][tid] = aco;
                    aro = Vcol[0][tid] * cos + Vcol[1][tid] * sin;
                    aco = Vcol[1][tid] * cos - Vcol[0][tid] * sin;
                    Vcol[0][tid] = aro;
                    Vcol[1][tid] = aco;

                    Aco[0][tid] = Acol[0][tid];
                    Vco[0][tid] = Vcol[0][tid];
                    tid = tid + blockDim.x;
                }
                __syncthreads();
            }
            else{
                tid = threadIdx.x;
                while (tid < row) {
                    aro = Acol[0][tid] * cos + Acol[1][tid] * sin;
                    aco = Acol[1][tid] * cos - Acol[0][tid] * sin;
                    Acol[0][tid] = aro;
                    Acol[1][tid] = aco;
                    aro = Vcol[0][tid] * cos + Vcol[1][tid] * sin;
                    aco = Vcol[1][tid] * cos - Vcol[0][tid] * sin;
                    Vcol[0][tid] = aro;
                    Vcol[1][tid] = aco;

                    Aco[1][tid] = Acol[1][tid];
                    Vco[1][tid] = Vcol[1][tid];
                    tid = tid + blockDim.x;
                }
                __syncthreads();
            }

            if (threadIdx.x == 0) {
                d1[pa[0]] = ele1 + ele * tan;
                d1[pa[1]] = ele2 - ele * tan;
            }
            __threadfence();
        }
		else if ( p_index == 0 ){
            tid = threadIdx.x;
            while (tid < row) {
                Aco[0][tid] = Acol[0][tid];
                Vco[0][tid] = Vcol[0][tid];
                tid = tid + blockDim.x;
            }
            __syncthreads();
		}

		__gpu_syncNew(Arrayin, i+1);
		//__gpu_sync(Arrayin);
	}

	if ( p_index == 0 ){
        tid = threadIdx.x;
        while (tid < row) {
            Aco[1][tid] = Acol[1][tid];
            Vco[1][tid] = Vcol[1][tid];
            tid = tid + blockDim.x;
        }
        __syncthreads();
    }
    else{
        tid = threadIdx.x;
        while (tid < row) {
            Aco[0][tid] = Acol[0][tid];
            Vco[0][tid] = Vcol[0][tid];
            tid = tid + blockDim.x;
        }
        __syncthreads();
    }
}

__device__ void updateMatrixA (int columnIndex[50], double *A, double *V, int row, int col, size_t pitch, double JacobiMatrix[50][50], int columnIndexJacobi[50][50], int numrowJacobiM[50], int numColumnJacobiM){

    int tid;
    double temp1, temp2;

    tid=threadIdx.x;
    while ( tid < row ){

        double *Arow=(double*)((char*)A+tid*pitch);
        double *Vrow=(double*)((char*)V+tid*pitch);

        for ( int j = 0; j < numColumnJacobiM; j ++ ){

            temp1 = 0.0;
            temp2 = 0.0;

            for ( int i = 0; i < numrowJacobiM[j]; i ++ ){
                temp1 += Arow[columnIndexJacobi[i][j]] * JacobiMatrix[i][j];
                temp2 += Vrow[columnIndexJacobi[i][j]] * JacobiMatrix[i][j];
            }

            Arow[columnIndex[j]] = temp1;
            Vrow[columnIndex[j]] = temp2;
        }

        tid+=blockDim.x;
    }
    __syncthreads();
}

__device__ void updateJacobiMatrix(double cosTheta, double sinTheta, int index1, int index2, double JacobiMatrix[50][50], int columnIndexJacobi[50][50], int numrowJacobiM[50], int numColumnJacobiM) {

    __shared__ double temp[50][2];
    __shared__ double indexRow[50];
    __shared__ int number;

    int tid;
    tid = threadIdx.x;
    while ( tid < numrowJacobiM[index1] ){
        temp[tid][0] = JacobiMatrix[tid][index1] * cosTheta;
        temp[tid][1] = JacobiMatrix[tid][index1] * sinTheta;

        indexRow[tid] = columnIndexJacobi[tid][index1];

        tid+=blockDim.x;
    }
    __syncthreads();

    if ( threadIdx.x == 0 ){

        number = numrowJacobiM[index1];

        for ( int j = 0; j < numrowJacobiM[index2]; j ++ ){

            int kk = 0;
            for ( int i = 0; i < numrowJacobiM[index1]; i ++ ){
                if ( columnIndexJacobi[i][index1] == columnIndexJacobi[j][index2] ){
                    kk = 1;
                    temp[i][0] += JacobiMatrix[j][index2] * sinTheta;
                    temp[i][1] -= JacobiMatrix[j][index2] * cosTheta;
                    break;
                }
            }

            if ( kk == 0 ){
                temp[number][0] =  1.0 * JacobiMatrix[j][index2] * sinTheta;
                temp[number][1] = -1.0 * JacobiMatrix[j][index2] * cosTheta;

                indexRow[number] = columnIndexJacobi[j][index2];

                number += 1;
            }
        }

        numrowJacobiM[index1] = number;
        numrowJacobiM[index2] = number;
    }
    __syncthreads();

    tid = threadIdx.x;
    while ( tid < number ){
        JacobiMatrix[tid][index1] = temp[tid][0];
        JacobiMatrix[tid][index2] = temp[tid][1];

        columnIndexJacobi[tid][index1] = indexRow[tid];
        columnIndexJacobi[tid][index2] = indexRow[tid];

        tid += blockDim.x;
    }
    __syncthreads();
}

__device__ void computeInnerProduct(double *cache, int index1, int index2, double *A, int row, int col, size_t pitch, double JacobiMatrix[50][50], int columnIndexJacobi[50][50], int numrowJacobiM[50], int numColumnJacobiM){

    cache[threadIdx.x]=0.0;

    double temp1, temp2;

    int tid=threadIdx.x;
    while( tid < row ){

        temp1 = 0.0;
        temp2 = 0.0;

        double *Arow=(double*)((char*)A+tid*pitch);

        for ( int i = 0; i < numrowJacobiM[index1]; i ++ ){
            temp1 += Arow[columnIndexJacobi[i][index1]] * JacobiMatrix[i][index1];
        }

        for ( int i = 0; i < numrowJacobiM[index2]; i ++ ){
            temp2 += Arow[columnIndexJacobi[i][index2]] * JacobiMatrix[i][index2];
        }

        cache[threadIdx.x] += temp1*temp2;
        tid+=blockDim.x;
    }
    __syncthreads();

    for( int i = 2; i <= blockDim.x; i=i*2 ){

        if ( threadIdx.x < blockDim.x/i ){
            cache[threadIdx.x]+=cache[threadIdx.x+(blockDim.x/i)];
        }
        __syncthreads();
    }
}

__global__ void gpusvdNew(double *A, double *V, int row, int col, int totalNumColumnGroup, int totalNumElemInGroup, size_t pitch,int *pass,int *count_i) {

    //int groupColumnIndex[2*gridDim.x][col/(2*gridDim.x)];
    //int totalNumColumnGroup = 2*gridDim.x;
    //int totalNumElemInGroup = col/(2*gridDim.x);

    __shared__ double ele, ele1, ele2;
    __shared__ double cache[256];

    __shared__ double JacobiMatrix[50][50];
    __shared__ int columnIndexJacobi[50][50];
    __shared__ int numrowJacobiM[50];
    __shared__ int numColumnJacobiM;
    int columnIndex[50];
    int localIndexGroup[2];
    int localIndexColumn[2];
    int tid, index1, index2;

    numColumnJacobiM = col/gridDim.x;

    for ( int i = 0; i < 2*gridDim.x-1; i ++ )
    {

        // get two local indexes of column groups
        localIndexGroup[0] = p_ab[0][blockIdx.x][i];
        localIndexGroup[1] = p_ab[1][blockIdx.x][i];

        // initialize Jacobi matrix
        tid = threadIdx.x;
        while ( tid < numColumnJacobiM ){
                JacobiMatrix[0][tid] = 1.0;
                numrowJacobiM[tid] = 1;
                tid += blockDim.x;
        }
        __syncthreads();

        tid = threadIdx.x;
        while ( tid < totalNumElemInGroup ){
            columnIndexJacobi[0][tid] = p1[groupColumnIndex[localIndexGroup[0]][tid]];
            columnIndexJacobi[0][tid+totalNumElemInGroup] = p1[groupColumnIndex[localIndexGroup[1]][tid]];

            columnIndex[tid] = p1[groupColumnIndex[localIndexGroup[0]][tid]];
            columnIndex[tid+totalNumElemInGroup] = p1[groupColumnIndex[localIndexGroup[1]][tid]];

            tid += blockDim.x;
        }
        __syncthreads();


        if ( i == 0 ){

            for ( int j = 0; j < 2*totalNumElemInGroup; j ++ ){

                // get the local index of Column 0
                localIndexColumn[0] = p1[blockIdx.x*(2*totalNumElemInGroup)+j];
                index1 =  totalNumElemInGroup * (j % 2) + j/2;

                for ( int k = j+1; k < 2*totalNumElemInGroup; k ++ ) {

                    // get the local index of Column 1
                    localIndexColumn[1] = p1[blockIdx.x*(2*totalNumElemInGroup)+k];
                    index2 = totalNumElemInGroup * (k % 2) + k/2;

                    // compute the inner product of two columns
                    computeInnerProduct(cache, index1, index2, A, row, col, pitch, JacobiMatrix, columnIndexJacobi, numrowJacobiM, numColumnJacobiM);

                    if(threadIdx.x==0){
                        ele=cache[0];
                    }
                    __syncthreads();

                    // compute cosTheta and sinTheta
                    double cosTheta = 0.0;
                    double sinTheta = 1.0;

                    if (fabs(ele)>THREASHOLD) {

                        if (threadIdx.x == 0) {
                            ele1 = d1[localIndexColumn[0]];
                            ele2 = d1[localIndexColumn[1]];
                            pass[blockIdx.x] = 0;
                        }
                        __syncthreads();

                        double tao = (ele1 - ele2) / (2 * ele);
                        double tan = sign(tao) / (fabs(tao) + sqrt(1 + pow(tao, 2)));
                        cosTheta = 1 / sqrt(1 + pow(tan, 2));
                        sinTheta = cosTheta * tan;

                        // update norms of two columns
                        if (threadIdx.x == 0) {
                            d1[localIndexColumn[0]] = ele1 + ele * tan;
                            d1[localIndexColumn[1]] = ele2 - ele * tan;
                        }
                        __threadfence();

                        // update Jacobi matrix
                        updateJacobiMatrix(cosTheta, sinTheta, index1, index2, JacobiMatrix, columnIndexJacobi, numrowJacobiM, numColumnJacobiM);
                    }
                }
            }
        }
        else {

            for ( int j = 0; j < totalNumElemInGroup; j ++ ){

                // get the local index of Column 0
                localIndexColumn[0] = columnIndex[j]; //p1[groupColumnIndex[localIndexGroup[0]][j]];
                index1 = j;

                for ( int k = 0; k < totalNumElemInGroup; k ++ ) {

                    // get the local index of Column 1
                    localIndexColumn[1] = columnIndex[totalNumElemInGroup+k]; // p1[groupColumnIndex[localIndexGroup[1]][k]];
                    index2 = totalNumElemInGroup + k;

                    // compute the inner product of two columns
                    computeInnerProduct(cache, index1, index2, A, row, col, pitch, JacobiMatrix, columnIndexJacobi, numrowJacobiM, numColumnJacobiM);

                    if(threadIdx.x==0){
                        ele=cache[0];
                    }
                    __syncthreads();

                    // compute cosTheta and sinTheta
                    double cosTheta = 0.0;
                    double sinTheta = 1.0;

                    if (fabs(ele)>THREASHOLD) {
                        if (threadIdx.x == 0) {
                            ele1 = d1[localIndexColumn[0]];
                            ele2 = d1[localIndexColumn[1]];
                            pass[blockIdx.x] = 0;
                        }
                        __syncthreads();

                        double tao = (ele1 - ele2) / (2 * ele);
                        double tan = sign(tao) / (fabs(tao) + sqrt(1 + pow(tao, 2)));
                        cosTheta = 1 / sqrt(1 + pow(tan, 2));
                        sinTheta = cosTheta * tan;

                        // update norms of two columns
                        if (threadIdx.x == 0) {
                            d1[localIndexColumn[0]] = ele1 + ele * tan;
                            d1[localIndexColumn[1]] = ele2 - ele * tan;
                        }
                        __threadfence();

                        // update Jacobi matrix
                        updateJacobiMatrix(cosTheta, sinTheta, index1, index2, JacobiMatrix, columnIndexJacobi, numrowJacobiM, numColumnJacobiM);
                    }
                }
            }
        }

        __syncthreads();

        // update Matrix A and V
        updateMatrixA (columnIndex, A, V, row, col, pitch, JacobiMatrix, columnIndexJacobi, numrowJacobiM, numColumnJacobiM);

        __gpu_syncNew(Arrayin, i+1);
    }
}


__global__ void getGroupColumnIndex ( int totalNumColumnGroup,  int totalNumElemInGroup ){

    int tid = threadIdx.x;

    while ( tid < totalNumColumnGroup ){
        for ( int i = 0; i < totalNumElemInGroup; i ++ ){
            groupColumnIndex[tid][i] = 2 * totalNumElemInGroup * (tid/2) + 2 * i + tid % 2;
        }
        tid += blockDim.x;
    }
    __syncthreads();

    /*
    if ( threadIdx.x == 0 ){

        for ( int j = 0; j < totalNumColumnGroup; j ++ ){
            printf("numIndex of group = %d\n", j);
            for ( int i = 0; i < totalNumElemInGroup; i ++ ){
                printf("%d \t", groupColumnIndex[j][i]) ;
            }
            printf("\n");
        }
    }
    __syncthreads();
    */
}

int main(int argc, char*argv[]){
	//int *dev_procSize;
	int width,height;
	FILE *fp;
    fp=fopen(argv[1],"r");//���ļ�
	fscanf(fp,"%d",&height);
	fscanf(fp,"%d",&width);
	double *host_A=(double*)malloc(width*height*sizeof(double));
	double *host_V=(double*)malloc(height*height*sizeof(double));
    for(int i=0;i<height;i++) //������
        for(int j=0;j<height;j++)
            fscanf(fp,"%lf",&host_A[i*width+j]);
    fclose(fp);//�ر��ļ�
	clock_t t1 = clock(); //��ʼ��ʱ��
	procSize = 125; //height/(COL*2);
    double *dev_A,*dev_A1,*dev_S,*dev_U,*dev_V,*dev_V1;
    for(int j=0;j<height;++j){
        host_V[j*height+j]=1.0;
    }
	double *host_U=(double*)malloc(width*width*sizeof(double));
	double *host_S=(double*)malloc(height*sizeof(double));
	size_t pitch;
	int *pass,*dev_pass,*count,*dev_count;
	cudaEvent_t start1;  
	cudaEventCreate(&start1);  
	cudaEvent_t stop1;  
	cudaEventCreate(&stop1);
	cudaHostAlloc((void**)&pass,height*sizeof(int),cudaHostAllocWriteCombined|cudaHostAllocMapped);
	cudaHostAlloc((void**)&count,sizeof(int),cudaHostAllocMapped);
	cudaHostGetDevicePointer(&dev_pass,pass,0);
	cudaHostGetDevicePointer(&dev_count,count,0);
    cudaMallocPitch(&dev_A,&pitch,width*sizeof(double),height);
    cudaMallocPitch(&dev_V,&pitch,height*sizeof(double),height);
	cudaMallocPitch(&dev_A1,&pitch,width*sizeof(double),height);
	cudaMallocPitch(&dev_V1,&pitch,height*sizeof(double),height);
	cudaMalloc((void**)&dev_U,width*width*sizeof(double));
	cudaMalloc((void**)&dev_S,height*sizeof(double));
    cudaMemcpy2D(dev_A,pitch,host_A,width*sizeof(double),width*sizeof(double),height,cudaMemcpyHostToDevice);
    cudaMemcpy2D(dev_V,pitch,host_V,height*sizeof(double),height*sizeof(double),height,cudaMemcpyHostToDevice);
	cudaEventRecord(start1, NULL);

    getRankNewNew<<<1,procSize>>>(2*procSize);
    int gridDim = procSize;
    int totalNumColumnGroup = 2*gridDim;
    int totalNumElemInGroup = width/(2*gridDim);

    getGroupColumnIndex<<<1,procSize>>>(totalNumColumnGroup, totalNumElemInGroup);  //[2*gridDim][col/(2*gridDim)];)

	int iteration= ITERATION;

    //while(iteration-- > 0)
    {
        //printf("iter = %d\n", iteration);
        count[0]=0;
        getD<<<height,256>>>(dev_A,height,width,pitch);
        cudaDeviceSynchronize();
        getrank<<<1,256>>>(height,width,dev_pass);
        cudaDeviceSynchronize();
        //getrankNew<<<1,width/2>>>(height,width,dev_pass);
        //cudaDeviceSynchronize();
        //getD<<<height,256>>>(dev_A,height,width,pitch);
        //cudaDeviceSynchronize();
        //getD<<<height,256>>>(dev_A,height,width,pitch);
        //cudaDeviceSynchronize();
        //gpusvd<<<procSize,256>>>(dev_A,dev_V,height,width,pitch,dev_pass,dev_count);
        gpusvdNew<<<procSize,256>>>(dev_A, dev_V, height, width, totalNumColumnGroup, totalNumElemInGroup, pitch,dev_pass, dev_count);
        cudaDeviceSynchronize();
        int allpass=0;
        for(int i=0;i<procSize;i++){
            allpass=allpass+pass[i];
        }

        if(allpass==procSize){
           //break;
        }
    }

	sort1<<<1,256>>>(dev_A,dev_V,dev_U,dev_S,dev_A1,dev_V1,height,width,pitch);
	cudaDeviceSynchronize();
	getU<<<1,256>>>(dev_A1,dev_S,dev_U,height,width,pitch);
	cudaDeviceSynchronize();
	printf("number of iteration=%d\n",ITERATION-iteration-1);
	cudaEventRecord(stop1, NULL);  
	cudaEventSynchronize(stop1);  
	float msecTotal1 = 0.0f;  
	cudaEventElapsedTime(&msecTotal1, start1, stop1); 
	printf("time = %f\n",msecTotal1);
	cudaError_t  error_check101 ;
			error_check101 = cudaGetLastError();
	        if( error_check101 != cudaSuccess ){
            printf("%s\n" , cudaGetErrorString( error_check101 ) );
         }
	cudaMemcpy2D(host_V,height*sizeof(double),dev_V1,pitch,height*sizeof(double),height,cudaMemcpyDeviceToHost);
	cudaMemcpy(host_S,dev_S,height*sizeof(double),cudaMemcpyDeviceToHost);
	cudaMemcpy(host_U,dev_U,width*width*sizeof(double),cudaMemcpyDeviceToHost);
	clock_t t2 = clock(); //������ʱ
	char outFileV[7]={'V','X','.','m','a','t','\0'};
    writeToFile(host_V,height,height,outFileV);
	char outFileU[7]={'U','X','.','m','a','t','\0'};
    writeToFile(host_U,width,width,outFileU);
	char outFileS[7]={'S','X','.','m','a','t','\0'};
    vectorToFile(host_S,height,outFileS);
    printf("cputime = %lf\n",(double)((t2-t1)*1000/CLOCKS_PER_SEC));
	cudaFree(dev_A);
	cudaFree(dev_V);
	cudaFree(dev_A1);
	cudaFree(dev_V1);
	cudaFree(dev_U);
	cudaFree(dev_S);
	cudaFree((void*)count);
	cudaFree((void*)pass);
	free(host_A);
	free(host_V);
	free(host_U);
	free(host_S);
    return 0;
}